# ------------------------------------------------------------------------------------------------
### Generate Figure 9: Map of LGBTQ-free zones for the cutoff sample
# ------------------------------------------------------------------------------------------------

# Load data
load("../data/processed/LGBTQ_election_data.Rda")

# Define color scheme and labels for the map
map_colors <- c("grey", "#D9534F", "white")
map_labels <- c(
  "No resolution < 50 km from the boundary", 
  "Anti-LGBTQ resolution < 50 km from the boundary", 
  "Beyond 50 km from the boundary"
)

# Load and process municipality boundaries
PL_gmina <- gisco_get_lau(country = "PL", year="2020") %>%
  mutate(id = paste0(substr(LAU_ID, 5, 6), substr(LAU_ID, 10, nchar(LAU_ID)))) %>%
  left_join(
    LGBTQ_election_data %>% 
      select(id, any_level_any_type) %>% 
      unique(), 
    by = "id"
  )

# Unite polygons within treated and control municipalities
united_treated <- PL_gmina %>% 
  filter(any_level_any_type == 1) %>%
  st_make_valid() %>%
  st_union()

united_control <- PL_gmina %>% 
  filter(any_level_any_type == 0) %>%
  st_make_valid() %>%
  st_union()

# Convert united polygons into MultiLineStrings and create policy adoption boundary
boundary_treated <- st_boundary(united_treated)
boundary_control <- st_boundary(united_control)
cutoff <- st_intersection(boundary_treated, boundary_control)

# Calculate municipal centroids and distances to the boundary
PL_gmina <- PL_gmina %>%
  mutate(
    centroid = st_centroid(`_ogr_geometry_`),
    dist = st_distance(centroid, cutoff),
    dist_km = as.numeric(dist / 1000)
  )

# Find municipalities within 50 km from the policy boundary
cutoff_gminas <- PL_gmina %>% 
  filter(dist_km < 50)

# Generate map
ggplot(PL_gmina) +
  geom_sf(aes(fill = "white")) +
  geom_sf(data = cutoff_gminas, aes(fill = any_level_any_type)) +
  geom_sf(data = cutoff, fill = NA, color = "black", size = 2.5) +
  scale_fill_manual(
    values = map_colors,
    labels = map_labels
  ) +
  theme_map() +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    text = element_text(size = 23)
  ) +
  guides(fill = guide_legend(ncol = 1))

# Save the plot
ggsave("../output/plots/figure9.pdf", width = 14, height = 10, dpi = 700)
